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A  Stability  Analysis  of  the  Perfectly  Matched  Layer  Method 

S.  Joe  Yakura,  David  Dietz,  Andy  Greenwood  and  Ernest  Baca 

*  Air  Force  Research  Laboratory,  Directed  Energy  Directorate, 

Kirtland  AFB,  New  Mexico  87117 


Abstract 

We  perform  a  detailed  stability  analysis  based  on  the  unsplit-field  uniaxial  perfectly  matched  layer  (PML) 
formulation .  Our  finding  shows  that  it  is  essential  to  have  transverse  field  gradients  present  at  all  times  to 
stabilize  PML  calculations .  In  the  absence  of  transverse  field  gradients,  the  PML  method  becomes  unstable 
with  the  axial  field  components  growing  linearly  in  time. 

L  INTRODUCTION 

Despite  the  successful  implementation  of  the  perfectly  matched  layer  (PML)  method  to  absorb  outgoing 
waves  at  the  artificial  boundaries  of  a  bounded  numerical  volume,  the  question  of  the  stability  of  the  PML 
method  remains  [1,2,3].  Abarbanel  and  Gottlieb  [1]  carried  out  a  detailed  stability  analysis  of  Berenger’s 
split-field  PML  formulation  [4],  and  they  concluded  that  the  split-field  PML  equations  are  not 
mathematically  strongly  well-posed.  Hence,  these  equations  result  in  unstable  field  components  that 
diverge  linearly  in  time. 

In  this  paper  we  present  a  stability  analysis  starting  with  the  unsplit-field  uniaxial  PML  formulation 
[5,6],  and  derive  a  stability  condition  for  the  simple,  nondispersive  PML  equations.  The  analysis  shows  that 
in  rare  instances  the  PML  method  results  in  an  unstable  condition.  However,  for  PML  parameter  values 
used  in  most  practical  applications,  the  PML  method  is  stable. 

II.  MAXWELL’S  EQUATIONS  INSIDE  PML 

For  a  plane  wave  propagating  in  the  arbitrary  x-y  direction  into  a  uniaxial  anisotropic  medium,  the  2-D 
PML  equations  in  the  frequency  (co)  domain  are  given  by  [5,6] 

V  xH(co;  x,  y)  =  jcoeo£RSPML(o);  x)E(co;  x,  y) ,  (II.  1 ) 

VxE((o;x,y)  =  -  jco)j.o|iRSPML  (ox  x)  H(co;  x,  y) ,  (H2) 

where  E(a>;x,y),  H(co;x,y),  eo,  £r,  Mo,  |Xr  and  SPML(co)  are  the  electric  field  vector,  the  magnetic  field  vector, 
the  free-space  permittivity,  the  relative  permittivity,  the  free-space  permeability,  the  relative  permeability 
and  the  uniaxial  anisotropic  PML  tensor,  respectively.  Elements  of  the  uniaxial  anisotropic  PML  tensor, 
SPML(co),  are  given  by 
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SPML(w)  = 


'  1 
Sx«o) 
0 
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0  0 

Sx(co)  0 

0  Sx(co) 


(H3) 


where  Sx(co)  is  an  arbitrarily  defined  co  and  x  dependent  function.  It  is  a  common  practice  in  the  FDTD 
community  to  choose  Sx(co)  in  the  form 
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with 
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(11*4) 


for  the  PML  matching  condition.  The  quantities  <rx  and  ax*  represent  electric  and  magnetic  conductivities 
arbitrarily  introduced  in  order  to  implement  the  PML  method. 

Since  Maxwell’s  equations  in  2-D  decompose  into  transverse  magnetic  (TMZ)  and  transverse  electric 
(TEZ)  waves,  these  waves  are  considered  separately  in  the  stability  analysis. 

II. A  TRANSVERSE  MAGNETIC  (TMZ)  WAVE 

For  TM2  waves,  Eqs.  (II.  1)  through  (II.4)  reduce  to  the  following  three  equations  for  three  field 
components  Hx(co;x,y),  Hy(co;x,y)  and  Ez(co;x,y): 


„  av 

[V  xH(co;  x,  y)]z  =  j(0£o£R  ( 1 4- - )EZ  (co;  x,  y) , 

jCO£o£R 

[V  xE(w,x,y)]x  = - Hx((0;x, y), 

(1  +  T^5-) 

jcoeogR 

[V  xE(^x,y)]  =  -jmpop.R(l+— )H  (to;x,y) . 

3  jcoeoER  3 


(n.5) 

(II.6) 

(n.7) 


Taking  the  inverse  Fourier  transform  of  the  above  equations  results  in  the  following  time-dependent  forms: 

(n.8) 


3Ez(t;x,y)  ox  1  3Hx(t;x,y)  1  3Hy(t;x,y) 
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0 1  £o£r  £o£r  d  y  £o£r  d  X 


8H,(l;x.y)  y)+ 

d  t  £o£r 


^V»(t;x,y)+-L®4^  =  0, 


^  £o£r  j 


jLlo|LiR  d  y 


3VYH(t;x,y)  aY  H 

a  VxH(t;x,y)-Hx(t;x,y)  =  0, 

d  t  £o£r 

(t:x,y)__!_8MiM>=0. 


d  t  &i£r 


JlIo(1r  3x 


(119) 

(11.10) 

(11.11) 


In  the  above,  the  first  order  time-dependent  equation  [Eq.  (II.  10)]  for  V„H(t;x,y)  is  introduced  to  handle  the 
delayed  time-response  of  Hx(t;x,y).  This  equation  follows  naturally  from  the  inverse  Fourier  transform  of 
Hx( oxx.y)  when  VxH(t;x,y)  is  defined  in  the  following  integral  form: 

t 

VxH(t;x,y)=  [Hx(t';x,y)  exp[-(-^-)(t-t')]  dt' .  (II. 1 2) 

J  £o£r 

— OO 

Casting  Eqs.  (II. 8)  through  (II.  1 1)  into  a  more  compact  form  results  in 
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(11.13) 


where  W™  =  (Ez(t;x,y),  Hx(t;x,y),  VxH(t;x,y),  Hy(t;x,y)  }T,  •  is  used  to  denote  matrix  multiplication,  and 
matrix  coefficients  A™  ,  B™  and  C™  are  given  by 
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(11.15-16) 


To  carryout  the  stability  analysis  of  Eq.  (11.13),  the  Laplace  transform  is  performed  in  the  time  domain 
and  the  Fourier  transform  in  the  spatial  domain.  The  Laplace  and  Fourier  transform  function  W™(s;kx,kv) 
of  W™(t;x,y)  is  defined  by 


S0+J« 


W™(t;x,y)— — — 7  Jexp(st)ds — — —  J  Jw™(s;kx,ky)exp(jkxx  +  jkyy)dkxdky  with  Re(s)>s0. 

(2jc)~  — oo 


Upon  performing  both  Laplace  and  Fourier  transforms,  Eq.  (11.13)  becomes 


(11.17) 


2™  (s; k x ,  k y )  •  W™  (s;  k x ,  k y )  =  W™  (0;  k x ,  k  ) , 


(11.18) 


where  U ™  (s;kx,ky)  is  the  characteristic  matrix  of  the  TMZ  wave  defined  by 


2™(s;kx,ky)  = 
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The  stability  of  the  system  is  characterized  by  investigating  the  determinant  (or  equivalently  the 
eigenvalues)  of  the  characteristic  matrix  D™  (s;kx,ky).  The  determinant  gives  the  following  quartic 
algebraic  equation: 
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^4+joi 
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(11.20) 


II.B  TRANSVERSE  ELECTRIC  (TEZ)  WAVE 

For  TEZ  waves,  Eqs.  (II.  1)  through  (II.4)  reduce  to  the  following  three  equations  for  three  field 
components  Hz(co;x,y),  Ex(co;x,y)  and  Ey(o);x,y): 


[V  xE(co;  X,  y)]z  =  -  jco|iio}XR  ( 1 +- — — )HZ  (co;  x,  y) , 

JCD£o£r 

[V  xH(co;  x,  y)]x  =  jt°£°£R  Ex(co;  x,  y) , 

(1+T^5-) 

JCQ£o£r 


(11.21) 


(11.22) 


[V  xH((o;  x,  y)]  =  jco£oeR  ( 1  +- — - — )E(co;  x,  y) . 

y  JCO£o£r  y 


(11.23) 


Following  the  same  steps  as  in  the  TMZ  wave  case  yields  the  following  equations  for  the  TEZ  wave  in  the 
Laplace  and  Fourier  domains: 


TC  A  Tp  At,  TP 

(s;kx,k  )•  W  (s;kx,ky)  =  W  (0;kx,ky), 


(11.24) 


where  WTE(s;kx,ky)  =  { Hz(s;kx,ky),  Ex(s;kx,ky),  VxE(s;kx,ky),  Ey(s;kx,ky)  }T,  and  Ote  (s;kx,ky)  is  the 
characteristic  matrix  of  the  TEZ  wave  defined  by 


£o£r  JLioJLiR 


Q  (s;kx,kv)  =  £o£r 


(IL25) 


Taking  the  determinant  of  the  characteristic  matrix,  Qr  (s;kx,ky),  gives  the  following  characteristic 
equation  which  is  exactly  the  same  as  in  the  TMZ  wave  case: 


2  (  ox  f  (kx)2  (  cx  f  (ky)2  A 

sz  s+ — —  +  —  - —  +  s+ — —  - =  0. 

I  £o£r  I  £o£R{io(iR  1  £o£r  I  £o£r(J.o(J.r 


(11.26) 


HI .  STABILITY  ANALYSIS  FOR  A  UNIAXIAL  PML  MEDIUM 


To  study  Eq.  (11.20)  [or  Eq.  (11.26)]  we  first  normalize  s,  kx  and  ky  by  setting 
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^  £o£r  ^ 


V  J 

Now,  Eq.  (11.20)  [or  Eq.  (II.26)]  results  in  the  following  form: 
S2  [(S  +  l)2  +  (Kx)2]  +  (S+l)2(K  )2  =  0. 


(III.4) 


From  this  expression  it  is  immediately  apparent  that  if  Ky  is  zero  then  two  of  the  four  roots  are  located  at 
S  =  0  in  the  complex  S-plane,  which  results  in  unconditionally  unstable  behavior  that  grows  linearly  with 
time.  These  two  roots  are  associated  with  the  axial  field  component  and  the  delayed  time-response  function 
of  the  axial  field  component.  The  other  two  roots  are  related  to  the  incoming  and  outgoing  damped 
transverse  waves  that  propagate  as  exp  [-  (gx/£0£r)  t  ±  jkxx]. 

On  the  other  hand,  if  both  Kx  and  Ky  are  real  numbers  then  Eq.  (III.4)  has  four  complex  roots  [i.e.  two 
sets  of  complex  conjugate  roots],  and  that  the  real  parts  of  these  roots  can  be  shown  to  be  all  negative  [see 
Appendix  A].  Thus,  all  eigenfunctions  associated  with  these  eigenvalues  are  well-behaved  and  stable. 

As  seen  in  Eq.  (III.4),  the  term  that  contributes  to  stabilizing  the  system  is  the  real  part  of  Ky. 
Physically,  this  means  that  the  transverse  field  gradients  (in  the  y-direetion  for  the  present  analysis) 
contribute  to  stabilize  axial  field  (in  the  x-direction  for  the  present  analysis)  components  as  TMZ  and  TEZ 
waves  propagate  into  a  uniaxial  anisotropic  PML  medium. 

Unfortunately,  the  actual  PML  system  is  not  typically  characterized  by  real  Kx  but  rather  by  complex  Kx 
because  of  the  evanescent  behavior  of  the  propagating  wave  into  the  uniaxial  PML  medium.  To  investigate 
the  effect  of  the  imaginary  part  of  Kx  on  the  stability  of  a  system,  we  solve  Eq.  (III.4)  directly  using 
MATHEMATICA  software  [7].  The  exact  expressions  for  the  four  complex  roots  are  shown  in 
Appendix  B.  Calculations  show  that  for  Im{  IKXI  }  »  Re{  IKXI  }  it  is  possible  for  the  real  parts  of  the  roots 
to  become  positive.  However,  in  usual  implementation  of  the  PML  method  Im{  IKXI  }  »  Re{  IKJ  }  is  not 
normally  satisfied;  thus,  it  is  unlikely  that  PML  calculations  become  unstable  for  practical  PML 
applications. 

IV  .  EXTENDING  THE  STABILITY  ANALYSIS  TO  CORNER  REGIONS 

At  a  corner  region  of  the  2D  PML  medium,  the  uniaxial  anisotropic  PML  tensor,  SPML(co),  has  to  be 
modified  to  include  contributions  in  both  x  and  y  directions  as  follows: 


SPML(oo) 
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where  Sy(co)  is  defined  in  the  same  way  as  Sx(to)  to  take  the  form 


(IV.  1) 


S  (oo)  =1  + - >_ 

jCOEoER 
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(IV  .2) 


Using  Eqs.  (IV.  1)  in  Eqs.  (II.  1)  and  (II.2),  and  following  the  same  steps  as  in  previous  sections  yields  the 
following  equations  for  the  TMZ  wave  in  the  Laplace  and  Fourier  domains  at  the  2D  PML  corner  region: 
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where  (  W™(s;kx,ky)  )corner  =  (Ez(s;kx,ky),  Hx(s;kx,ky),  VxH(s;kx,ky),  Hy(s;kx,ky),  VyH(s;kx,ky)  }T,  and 
(£2™  (s;kx,ky)  )comer  is  the  characteristic  matrix  of  the  TMZ  defined  by 
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(IV.4) 


Taking  the  determinant  of  the  characteristic  matrix,  (£2™  (s;kx,ky)  )  c01-nen  gives  the  following  characteristic 
equation: 
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(IV  .5) 


As  seen  in  the  above  equation  one  root  is  located  at  s  =  0,  which  gives  a  stable  solution  in  the  time 
domain,  and  the  other  four  roots  can  be  obtained  by  setting  the  expression  inside  the  square  bracket  to  zero. 
For  real  values  of  kx  and  ky,  a  procedure  similar  to  that  in  Appendix  A  shows  the  real  parts  of  the  four 
complex  roots  are  always  negative,  implying  stable  solutions  in  the  time  domain.  For  arbitrary  complex 
values  of  kx  and  ky,  the  real  parts  of  the  four  complex  roots  have  to  be  investigated  numerically  from  the 
exact  expressions  shown  in  Appendix  B. 

For  the  special  case  of  ox  =  c?y  the  equation  formed  by  setting  the  expression  inside  the  square  bracket  to 
zero  can  be  solved  exactly  to  obtain  an  analytical  expression  for  the  stability  condition;  in  this  case  the 
square  bracket  term  reduces  to 


\2 


s+- 


£o£r 


S  +  - 


^2  +  (kx)2+(kx)2 


£o£r 


£o£r]jLo}Ir 


=  0. 


(IV.6) 


Solving  Eq.  (IV.6)  results  in  the  double  root  s  =  -  (gx/£o£r)  ,  which  give  stable  solutions  in  the  time 

domain,  and  the  other  two  roots  s  =  -  (cx  /  £o£r)  ±  j  {^/l/ (£o£r|1o|lir)  k  }  where  (k)2  =  (kx)2  +  (ky)2.  Expressing 
k  in  terms  of  its  real  and  imaginary  parts  as  k=kR+jk!,  the  two  roots 
s  =  -  (ax  /£o£r)±  j{yl/(£o£Rp.ojiR)  (  kR  +  jk1 )  }  can  be  expressed  as 


£o£r 


±j 


[(kR)2+(k‘)2] 
£o£r|1o|Xr 


exp 


j  tan' 


/J 


(kK  *  0). 


(IV  .7) 
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In  the  above  expression,  one  of  the  two  roots  gives  a  positive  real  value  if  the  following  condition  is 
satisfied: 


Hence,  the  PML  system  becomes  unstable  if  the  above  condition  is  met  for  the  case  ox=  cry. 

The  stability  analysis  of  TEZ  waves  in  corner  regions  results  in  the  same  stability  condition  as  for  TMZ 
waves. 


V.  CONCLUSIONS 

Starting  with  unsplit-field  uniaxial  PML  formulation  in  the  frequency  domain.  Maxwell’s  equations  are 
cast  into  a  set  of  first  order  differential  equations  in  time.  Then,  using  the  Laplace  and  Fourier  transforms, 
the  characteristic  equation  of  a  system  is  obtained  and  investigated  for  its  dynamic  stability. 

From  stability  analysis,  we  find  that  it  is  essential  for  the  transverse  field  gradients  to  be  present  at  all 
times  in  order  to  stabilize  PML  calculations.  In  fact,  in  the  absence  of  transverse  field  gradients  the  PML 
method  becomes  unstable  with  the  axial  field  components  growing  linearly  in  time. 
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oendix  A 


A  proof  to  show  that  the  real  parts  of  the  roots  of  the  polynomial,  S2  [  (1  +  S  )2  +  a]  +  ( 1  +  S  )2  b  =  0  with 
positive  real  coefficients  a  and  b,  are  all  negative. 


Theorem: 

Consider  the  equation 


S2[(S  +  I)2+a]  +  (S  +  1)2b  =  0,  SeC 


(A.l) 
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where  a,  b  >  0.  Then 

(i)  Eq.  (A.  1 )  has  no  real  solutions; 

(ii)  If  Sj  =  oti  +  jPi,  i  =1, 2,  3, 4  (Pi ^  0)  are  the  roots  of  Eq.  (A.l)  then  oti  <  0,  i  =  1,  2,  3,  4. 

Proof: 

(i)  Rewrite  Eq.  (A.  1 )  as 

S2[(S+l)2  +a]  =  -(S+l)2b ,  (A.2) 

and  note  that,  since  a,  b  >  0,  if  S  G  R  \{-l,0}  then  the  LHS  of  Eq.  (A.2)  >  0  and  the  RHS  of  Eq.  (A.2)  <  0, 
a  contradiction;  while  if  S  =  -1  then  the  LHS  of  Eq.  (A.2)  >  0  and  the  RHS  of  Eq.  (A.2)  =  0  and  if  S  =  0 
then  the  LHS  of  Eq.  (A.2)  =  0  and  the  RHS  of  Eq.  (A.2)  <  0,  also  contradictions.  Thus,  if  S  satisfies 

Eq.  (A.l)  then  S  $  R.  Further,  all  four  solutions  are  of  the  form  S  =  a  +  j(J,  a,  p  €  R,  p  *  0. 

(ii)  If  S  =  a  +  jp  is  any  solution  of  Eq.  (A.l)  then 

(a  + jP)2[(a+ jP  +  l)2  +  a]  +  (a+ jP  +  l)2b  =  0 .  (A.3) 

Expanding  and  equating  real  and  imaginary  parts  of  Eq.  (A.3)  to  zero  gives 

(a2  -P2)2  -4a2p2  +  2[a  (a2  -p2)-2ap2]  +  (l  +  a  +  b)(a2  ~p2)  +  2ab  +  b  =  0  (A.4) 

and 

P[2ot(a2  -p2)  +  2a2  +(a2  -p2)  +  (l  +  a  +  b)a  +  b]  =  0 .  (A.5) 

In  Eq.  (A.5)  the  term  in  brackets  must  be  equal  to  zero  since,  otherwise,  p  =  0,  which  is  not  possible  since 
S  £  R.  Rewritten  the  term  results  in 

-(2a+l)p2  +  2a3  +3a2  +  (l  +  a  +  b)a  +  b  =  0 .  ■  (A.6) 

If  a  =  -Vi  then  we  are  done  (since  then  a  <  0).  Otherwise,  if  a  &  -Vi  then  Eq.  (A.6)  gives 


32  = — - — [2a3  +3a2  +(l  +  a  +  b)a  +  b] 

2a  +  1 

(A.7) 

P”  —  a-  —  [2a“  +  (l  +  a +b)a+ b] . 

2a+l  v  ’ 

(A.8) 

Substituting  Eqs.  (A.7)  and  (A. 8)  into  Eq.  (A.4)  and  simplifying  leads  to 

16cx6  +48  a5  +8(7  +  a +  b)  a4  +16(2  + a +  b)  a3 

+  [(l  +  a  +  b)2  +  8(l  +  a  +  b)]a2  +  (l  +  a  +  b)2a  +  ab  =  0.  (A.9) 

Since  a,  b  >  0  then  all  coefficients  in  Eq.  (A.9)  are  >  0;  thus,  by  Descartes's  Rule  of  Signs  [8],  Eq.  (A.9)  has 
no  positive  roots.  Further,  zero  is  not  a  root  of  Eq.  (A.9).  Hence,  all  real  roots  of  Eq.  (A.9)  are  <  0.  Finally, 
oti,  i  =  1, 2,  3,  4  must  be  among  the  solutions  of  Eq.  (A.9)  so  otj  <  0,  i  =  1, 2,  3,  4. 

Appendix  B 

Four  complex  roots  of  the  polynomial,  S4+aS3+bS2  +  cS  +  d  =0,  with  complex  coefficients  a,  b,  c  and  d 
are  given  by  [7] 
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